Sample Statistics
# Income (mean of both partner's report)
merge_income <- function(income1, income2) {
merged_income <- numeric(length(income1))
# Loop through each pair of incomes
for (i in seq_along(income1)) {
# Handle NA
if (is.na(income1[i])) {
merged_income[i] <- income2[i]
}
else if (is.na(income2[i])) {
merged_income[i] <- income1[i]
}
# if both are informative, take mean and round
else if (income1[i] %in% 1:6 && income2[i] %in% 1:6) {
merged_income[i] <- round((income1[i] + income2[i]) / 2)
}
# if one is informative and the other not, use the informative one.
else if (income1[i] %in% 1:6) {
merged_income[i] <- income1[i]
}
else if (income2[i] %in% 1:6) {
merged_income[i] <- income2[i]
}
# Now we only have cases, where both are either 70 or 99. We simply report "undisclosed".
else {
merged_income[i] <- 99
}
}
# Convert to factor
merged_income <- factor(
merged_income, levels = c(1,2,3,4,5,6,70,99),
labels = c(
"up to CHF 2'000.-",
"CHF 2'001.- to CHF 4'000.-",
"CHF 4'001.- to CHF 6'000.-",
"CHF 6'001.- to CHF 8'000.-",
"CHF 8'001.- to CHF 10'000.-",
"above CHF 10'000.-",
"I don't know",
"Undisclosed"
)
)
return(merged_income)
}
df_sample_report <- df_full %>%
group_by(coupleID) %>%
arrange(userID) %>%
# Computing couple level variables
mutate(
Household_Income = merge_income(first(pre_income_1), last(pre_income_1)),
reldur = pre_rel_duration_m / 12 + pre_rel_duration_y,
Relationship_duration = mean(reldur, na.rm = TRUE),
habdur = pre_hab_duration_m / 12 + pre_hab_duration_y,
Cohabiting_duration = mean(habdur, na.rm = TRUE),
Marital_status = factor(
case_when(
all(pre_mat_stat == 1) ~ "Married",
any(pre_mat_stat == 1) ~ "One Partner Married",
TRUE ~ "Not Married"
)
),
Have_children = factor(
(first(pre_child_option) + last(pre_child_option)) > 0,
levels = c(FALSE, TRUE),
labels = c('Have Children', 'No Children')),
Gender = factor(
gender,
levels = c(1,2,3),
labels = c('Male','Female', 'Other')),
Couple_type = as.factor(
case_when(
first(Gender) == last(Gender) & first(Gender) == 'Male' ~ 'Same-Sex Couple (Male)',
first(Gender) == last(Gender) & first(Gender) == 'Female' ~ 'Same-Sex Couple (Female)',
TRUE ~ 'Mixed-sex Couple'
)
)
) %>%
ungroup() %>%
# Individual level variables
mutate(
Age = pre_age,
Handedness = factor(
pre_handedness,
levels = c(0, 1, 2),
c('Right','Left', 'Ambidextrous')),
Highest_Education = factor(
pre_education,
levels = c(1,2,3,4,5,6,7),
labels = c(
"(still) no school diploma",
"compulsory education (9 years)",
"vocational training (apprenticeship)",
"Matura (university entrance qualification)",
"Bachelor's degree",
"Master's degree",
"Doctorate degree"
)
),
BMI = pre_weight / ((pre_height / 100)^2) # to meters
) %>%
select(c(Relationship_duration, Cohabiting_duration, Couple_type, Household_Income,
Marital_status, Have_children,
Gender, Age, Handedness, Highest_Education, BMI))
sample_table <- report_measures(df_sample_report, ICC = F)
sample_table$n_Obs <- as.numeric(sample_table$n_Obs) / 55
rownames(sample_table) <- NULL
n_couple_vars <- 17
sample_table$n_Obs[1:n_couple_vars] <- sample_table$n_Obs[1:n_couple_vars] / 2
packing_sample <- list(
"Couple level variables (38 couples)" =
c(1, n_couple_vars),
"Individual level variables (76 individuals)"
= c(n_couple_vars+1, nrow(sample_table))
)
df_sample_summary <- print_df(
sample_table,
rows_to_pack = packing_sample
)
export_xlsx(df_sample_summary,
file.path('Output', 'SampleDescription.xlsx'),
merge_option = 'both',
rows_to_pack = packing_sample,
colwidths = c(20,35,7,7,7,7,7,10)
)
##
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
##
## guess_encoding
|
Variable
|
Level
|
n_Obs
|
percentage_Obs
|
Missing
|
Mean
|
SD
|
Range
|
|
Couple level variables (38 couples)
|
|
Relationship_duration
|
NA
|
38
|
NA
|
0%
|
9.23
|
9.03
|
0.58-36.00
|
|
Cohabiting_duration
|
NA
|
38
|
NA
|
0%
|
7.53
|
9.14
|
0.25-33.00
|
|
Couple_type
|
Same-Sex Couple (Male)
|
1
|
3%
|
NA
|
NA
|
NA
|
NA
|
|
Couple_type
|
Same-Sex Couple (Female)
|
1
|
3%
|
NA
|
NA
|
NA
|
NA
|
|
Couple_type
|
Mixed-sex Couple
|
36
|
95%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
I don’t know
|
0
|
0%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
up to CHF 2’000.-
|
2
|
5%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
CHF 2’001.- to CHF 4’000.-
|
3
|
8%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
CHF 6’001.- to CHF 8’000.-
|
3
|
8%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
CHF 4’001.- to CHF 6’000.-
|
5
|
13%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
Undisclosed
|
5
|
13%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
CHF 8’001.- to CHF 10’000.-
|
8
|
21%
|
NA
|
NA
|
NA
|
NA
|
|
Household_Income
|
above CHF 10’000.-
|
12
|
32%
|
NA
|
NA
|
NA
|
NA
|
|
Marital_status
|
Married
|
13
|
34%
|
NA
|
NA
|
NA
|
NA
|
|
Marital_status
|
Not Married
|
25
|
66%
|
NA
|
NA
|
NA
|
NA
|
|
Have_children
|
No Children
|
11
|
29%
|
NA
|
NA
|
NA
|
NA
|
|
Have_children
|
Have Children
|
27
|
71%
|
NA
|
NA
|
NA
|
NA
|
|
Individual level variables (76 individuals)
|
|
Gender
|
Other
|
0
|
0%
|
NA
|
NA
|
NA
|
NA
|
|
Gender
|
Male
|
38
|
50%
|
NA
|
NA
|
NA
|
NA
|
|
Gender
|
Female
|
38
|
50%
|
NA
|
NA
|
NA
|
NA
|
|
Age
|
NA
|
76
|
NA
|
0%
|
34.01
|
10.96
|
19.00-60.00
|
|
Handedness
|
Ambidextrous
|
1
|
1%
|
NA
|
NA
|
NA
|
NA
|
|
Handedness
|
Left
|
11
|
14%
|
NA
|
NA
|
NA
|
NA
|
|
Handedness
|
Right
|
64
|
84%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
(still) no school diploma
|
0
|
0%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
compulsory education (9 years)
|
0
|
0%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
Doctorate degree
|
0
|
0%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
vocational training (apprenticeship)
|
8
|
11%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
Master’s degree
|
21
|
28%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
Bachelor’s degree
|
23
|
30%
|
NA
|
NA
|
NA
|
NA
|
|
Highest_Education
|
Matura (university entrance qualification)
|
24
|
32%
|
NA
|
NA
|
NA
|
NA
|
|
BMI
|
NA
|
76
|
NA
|
0%
|
24.94
|
4.08
|
16.37-33.95
|
Measures
Focal Variables
main_constructs <- c("persuasion", "pressure","pushing",
"pa_sub", "pa_obj", "aff", "reactance"
)
main_descriptives <- report_measures(
data = df_full,
measures = main_constructs,
ICC = TRUE,
cluster_var = df_full$userID)
openxlsx::write.xlsx(
main_descriptives,
file.path('Output', 'DescriptivesMain.xlsx')
)
print_df(main_descriptives)
|
Variable
|
n_Obs
|
Missing
|
Mean
|
SD
|
Range
|
ICC
|
|
persuasion
|
4180
|
6%
|
0.42
|
1.08
|
0.00- 5.00
|
0.23
|
|
pressure
|
4180
|
6%
|
0.12
|
0.62
|
0.00- 5.00
|
0.58
|
|
pushing
|
4180
|
6%
|
0.16
|
0.66
|
0.00- 5.00
|
0.11
|
|
pa_sub
|
4180
|
6%
|
30.40
|
54.78
|
0.00-720.00
|
0.16
|
|
pa_obj
|
4180
|
11%
|
144.41
|
117.81
|
5.75-971.25
|
0.18
|
|
aff
|
4180
|
6%
|
4.83
|
1.14
|
1.00- 6.00
|
0.41
|
|
reactance
|
4180
|
81%
|
0.79
|
1.32
|
0.00- 5.00
|
0.47
|
All Variables
all_constructs <- c(
main_constructs,
"day",
"weartime",
"isWeekend",
"plan",
"studyGroup",
"support",
"got_JITAI",
"skilled_support"
)
all_descriptives <- report_measures(df_full, all_constructs, ICC = F)
openxlsx::write.xlsx(
all_descriptives,
file.path('Output', 'DescriptivesAll.xlsx')
)
print_df(all_descriptives)
|
Variable
|
Level
|
n_Obs
|
percentage_Obs
|
Missing
|
Mean
|
SD
|
Range
|
|
persuasion
|
NA
|
4180
|
NA
|
6%
|
0.42
|
1.08
|
0.00- 5.00
|
|
pressure
|
NA
|
4180
|
NA
|
6%
|
0.12
|
0.62
|
0.00- 5.00
|
|
pushing
|
NA
|
4180
|
NA
|
6%
|
0.16
|
0.66
|
0.00- 5.00
|
|
pa_sub
|
NA
|
4180
|
NA
|
6%
|
30.40
|
54.78
|
0.00- 720.00
|
|
pa_obj
|
NA
|
4180
|
NA
|
11%
|
144.41
|
117.81
|
5.75- 971.25
|
|
aff
|
NA
|
4180
|
NA
|
6%
|
4.83
|
1.14
|
1.00- 6.00
|
|
reactance
|
NA
|
4180
|
NA
|
81%
|
0.79
|
1.32
|
0.00- 5.00
|
|
day
|
NA
|
4180
|
NA
|
0%
|
0.50
|
0.29
|
0.00- 1.00
|
|
weartime
|
NA
|
4180
|
NA
|
1%
|
1057.42
|
384.29
|
0.00-1440.00
|
|
isWeekend
|
Weekend
|
1216
|
29%
|
NA
|
NA
|
NA
|
NA
|
|
isWeekend
|
Weekday
|
2964
|
71%
|
NA
|
NA
|
NA
|
NA
|
|
plan
|
missing
|
238
|
6%
|
NA
|
NA
|
NA
|
NA
|
|
plan
|
Plan
|
1860
|
44%
|
NA
|
NA
|
NA
|
NA
|
|
plan
|
No plan
|
2082
|
50%
|
NA
|
NA
|
NA
|
NA
|
|
studyGroup
|
last 3 weeks interventions
|
1320
|
32%
|
NA
|
NA
|
NA
|
NA
|
|
studyGroup
|
Allways inerventions
|
1430
|
34%
|
NA
|
NA
|
NA
|
NA
|
|
studyGroup
|
First 3 weeks interventions
|
1430
|
34%
|
NA
|
NA
|
NA
|
NA
|
|
support
|
NA
|
4180
|
NA
|
6%
|
0.91
|
1.52
|
0.00- 5.00
|
|
got_JITAI
|
JITAI received
|
585
|
14%
|
NA
|
NA
|
NA
|
NA
|
|
got_JITAI
|
No JITAI
|
3595
|
86%
|
NA
|
NA
|
NA
|
NA
|
|
skilled_support
|
Days before Intervention
|
1036
|
25%
|
NA
|
NA
|
NA
|
NA
|
|
skilled_support
|
Days after Intervention
|
3144
|
75%
|
NA
|
NA
|
NA
|
NA
|
Correlations
cors <- wbCorr(df_full[,c(main_constructs)], df_full$coupleID, method = 'spearman')
main_cors <- summary(cors, 'wb')$merged_wb
openxlsx::write.xlsx(
main_cors,
file.path('Output', 'Correlations.xlsx')
)
print_df(main_cors, width = '7em')
|
|
persuasion
|
pressure
|
pushing
|
pa_sub
|
pa_obj
|
aff
|
reactance
|
|
persuasion
|
[0.20]
|
0.21***
|
0.40***
|
0.17***
|
0.07***
|
0.00
|
-0.05
|
|
pressure
|
0.30
|
[0.55]
|
0.28***
|
-0.03
|
-0.04*
|
-0.01
|
0.26***
|
|
pushing
|
0.63***
|
0.40*
|
[0.08]
|
0.14***
|
0.05**
|
0.01
|
0.07*
|
|
pa_sub
|
0.27
|
-0.10
|
0.24
|
[0.15]
|
0.31***
|
0.20***
|
-0.04
|
|
pa_obj
|
0.13
|
-0.08
|
0.27
|
0.51***
|
[0.13]
|
0.19***
|
0.06
|
|
aff
|
0.28
|
-0.07
|
0.29
|
0.52***
|
0.22
|
[0.28]
|
-0.01
|
|
reactance
|
0.18
|
0.23
|
0.11
|
0.06
|
0.38*
|
-0.12
|
[0.42]
|
Within-person correlations are above the diagonal and between-person
correlations are below the diagonal.
On the diagonal are intraclass correlations (ICCs)
Modelling
# For indistinguishable Dyads
model_rows_fixed <- c(
'Intercept',
# '-- WITHIN PERSON MAIN EFFECTS --',
'persuasion_self_cw',
'persuasion_partner_cw',
'pressure_self_cw',
'pressure_partner_cw',
'pushing_self_cw',
'pushing_partner_cw',
'day',
'weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'persuasion_self_cb',
'persuasion_partner_cb',
'pressure_self_cb',
'pressure_partner_cb',
'pushing_self_cb',
'pushing_partner_cb',
'weartime_self_cb'
)
model_rows_fixed_ordinal <- c(
model_rows_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rows_fixed[2:length(model_rows_fixed)]
)
model_rows_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(persuasion_self_cw)',
'sd(persuasion_partner_cw)',
'sd(pressure_self_cw)',
'sd(pressure_partner_cw)',
'sd(pushing_self_cw)',
'sd(pushing_partner_cw)',
# '-- CORRELATION STRUCTURE -- ',
'ar[1]',
'nu',
'shape',
'sderr',
'sigma'
)
model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
'Intercept',
# '-- WITHIN PERSON MAIN EFFECTS --',
'Daily received persuasion target -> target',
'Daily received persuasion target -> agent',
'Daily received pressure target -> target',
'Daily received pressure target -> agent',
'Daily received pushing target -> target',
'Daily received pushing target -> agent',
'Day',
'Daily weartime',
# '-- BETWEEN PERSON MAIN EFFECTS',
'Mean received persuasion target -> target',
'Mean received persuasion target -> agent',
'Mean received pressure target -> target',
'Mean received pressure target -> agent',
'Mean received pushing target -> target',
'Mean received pushing target -> agent',
'Mean weartime'
)
model_rownames_fixed_ordinal <- c(
model_rownames_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rownames_fixed[2:length(model_rownames_fixed)]
)
model_rownames_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(Daily received persuasion target -> target)',
'sd(Daily received persuasion target -> agent)',
'sd(Daily received pressure target -> target)',
'sd(Daily received pressure target -> agent)',
'sd(Daily received pushing target -> target)',
'sd(Daily received pushing target -> agent)',
# '-- CORRELATION STRUCTURE -- ',
'ar[1]',
'nu',
'shape',
'sderr',
'sigma'
)
model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
"Intercept",
# "-- WITHIN PERSON MAIN EFFECTS --",
"Daily persuasion experienced",
"Daily persuasion utilized (partner's view)", # OR partner received
"Daily pressure experienced",
"Daily pressure utilized (partner's view)",
"Daily pushing experienced",
"Daily pushing utilized (partner's view)",
"Day",
"Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Mean persuasion experienced",
"Mean persuasion utilized (partner's view)",
"Mean pressure experienced",
"Mean pressure utilized (partner's view)",
"Mean pushing experienced",
"Mean pushing utilized (partner's view)",
"Mean weartime"
)
model_rownames_fixed_ordinal <- c(
model_rownames_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rownames_fixed[2:length(model_rownames_fixed)]
)
model_rownames_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
"sd(Daily persuasion experienced)",
"sd(Daily persuasion utilized (partner's view))", # OR partner received
"sd(Daily pressure experienced)",
"sd(Daily pressure utilized (partner's view))",
"sd(Daily pushing experienced)",
"sd(Daily pushing utilized (partner's view))",
# '-- CORRELATION STRUCTURE -- ',
'ar[1]',
'nu',
'shape',
'sderr',
'sigma'
)
model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
"Within-Person Effects" = c(2,9),
"Between-Person Effects" = c(10,16),
"Random Effects" = c(17, 23),
"Additional Parameters" = c(24,28)
)
rows_to_pack_ordinal <- list(
"Intercepts" = c(1,6),
"Within-Person Effects" = c(2+5,9+5),
"Between-Person Effects" = c(10+5,16+5),
"Random Effects" = c(17+5, 23+5),
"Additional Parameters" = c(24+5,28+6)
)
Subjective MVPA
range(df_double$pa_sub, na.rm = T)
## [1] 0 720
hist(df_double$pa_sub, breaks = 100)

Modelling using the gaussian family fails. Due to the many zeros,
transformations won’t help estimating the models. We employ the negative
binomial family.
formula <- bf(
pa_sub ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw | coupleID),
autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b"),
brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
brms::set_prior("cauchy(0, 20)", class = "shape"),
brms::set_prior("cauchy(0, 10)", class='sderr')
#brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::negbinomial(),
#control = list(adapt_delta = 0.99),
iter = 5000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", "pa_sub")
)
plot(pa_sub, ask = FALSE)









plot(pp_check(pa_sub, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

pp_check_transformed(pa_sub, transform = log1p)

##
## Computed from 12000 by 3736 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -12060.6 177.6
## p_loo 30.2 2.8
## looic 24121.2 355.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.6]).
##
## Pareto k diagnostic values:
## Count Pct.
## (-Inf, 0.7] (good) 3732 99.9%
## (0.7, 1] (bad) 3 0.1%
## (1, Inf) (very bad) 1 0.0%
## Min. ESS
## (-Inf, 0.7] 662
## (0.7, 1] <NA>
## (1, Inf) <NA>
## See help('pareto-k-diagnostic') for details.
summarize_brms(
pa_sub,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)
|
|
IRR
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
Intercept
|
27.35*
|
21.03
|
35.95
|
1.001
|
3902.26
|
6229.26
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
1.20*
|
1.07
|
1.35
|
1.000
|
10192.93
|
7278.69
|
|
Daily persuasion utilized (partner’s view)
|
1.19*
|
1.06
|
1.34
|
1.000
|
12159.42
|
8804.61
|
|
Daily pressure experienced
|
0.94
|
0.71
|
1.30
|
1.000
|
9685.08
|
7640.07
|
|
Daily pressure utilized (partner’s view)
|
1.16
|
0.88
|
1.59
|
1.000
|
11249.31
|
7539.42
|
|
Daily pushing experienced
|
1.13
|
0.92
|
1.41
|
1.000
|
8816.39
|
8383.16
|
|
Daily pushing utilized (partner’s view)
|
1.13
|
0.95
|
1.36
|
1.000
|
12547.79
|
9213.16
|
|
Day
|
0.79
|
0.58
|
1.10
|
1.000
|
12431.32
|
9455.82
|
|
Daily weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
1.61
|
0.83
|
3.13
|
1.002
|
3172.48
|
5991.41
|
|
Mean persuasion utilized (partner’s view)
|
1.19
|
0.62
|
2.28
|
1.002
|
3128.46
|
5435.36
|
|
Mean pressure experienced
|
0.50
|
0.22
|
1.11
|
1.001
|
4626.65
|
6921.35
|
|
Mean pressure utilized (partner’s view)
|
0.43*
|
0.19
|
0.96
|
1.003
|
4346.50
|
7133.59
|
|
Mean pushing experienced
|
1.71
|
0.64
|
4.56
|
1.000
|
4336.92
|
5727.04
|
|
Mean pushing utilized (partner’s view)
|
2.03
|
0.76
|
5.67
|
1.001
|
4271.64
|
6953.82
|
|
Mean weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Random Effects
|
|
sd(Intercept)
|
0.61
|
0.45
|
0.83
|
1.00
|
3236.86
|
5929.19
|
|
sd(Daily persuasion experienced)
|
0.09
|
0.00
|
0.24
|
1.00
|
4277.64
|
4971.03
|
|
sd(Daily persuasion utilized (partner’s view))
|
0.08
|
0.00
|
0.21
|
1.00
|
5008.36
|
5206.23
|
|
sd(Daily pressure experienced)
|
0.17
|
0.01
|
0.53
|
1.00
|
6480.10
|
5580.93
|
|
sd(Daily pressure utilized (partner’s view))
|
0.16
|
0.01
|
0.49
|
1.00
|
7007.22
|
4748.59
|
|
sd(Daily pushing experienced)
|
0.28
|
0.02
|
0.58
|
1.00
|
3266.42
|
3678.71
|
|
sd(Daily pushing utilized (partner’s view))
|
0.12
|
0.00
|
0.32
|
1.00
|
4609.50
|
3727.29
|
|
Additional Parameters
|
|
ar[1]
|
0.03
|
-0.94
|
0.94
|
1.00
|
10005.25
|
6639.54
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
0.14
|
0.13
|
0.14
|
1.00
|
13938.19
|
8436.71
|
|
sderr
|
0.05
|
0.00
|
0.13
|
1.00
|
6131.05
|
5072.35
|
|
sigma
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
Device Based MVPA
range(df_double$pa_obj, na.rm = T)
## [1] 5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)
hist(df_double$pa_obj_log, breaks = 50)

We tried negative binomial here as well for consistency, but the
model did not converge. Poisson also did not work. As we have no zeros
in this distribution, we log transform.
formula <- bf(
pa_obj_log ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
pushing_self_cb + pushing_partner_cb +
day + weartime_self_cw + weartime_self_cb +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw | coupleID),
autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b"),
brms::set_prior("normal(0, 50)", class = "Intercept", lb = 0),
brms::set_prior("normal(0, 10)", class = "sd", group = "coupleID", lb = 0),
brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.99),
iter = 5000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", "pa_obj_log")
)
plot(pa_obj_log, ask = FALSE)










plot(pp_check(pa_obj_log, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(pa_obj_log))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(pa_obj_log, transform = log1p)
loo(pa_obj_log)
##
## Computed from 12000 by 3337 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -2813.7 56.5
## p_loo 91.6 4.5
## looic 5627.5 113.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 1.9]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
summarize_brms(
pa_obj_log,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)
|
|
exp(Est.)
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
Intercept
|
118.09*
|
106.05
|
131.83
|
1.002
|
4338.46
|
6795.25
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
1.03
|
1.00
|
1.06
|
1.001
|
10995.72
|
9537.78
|
|
Daily persuasion utilized (partner’s view)
|
1.02
|
0.99
|
1.05
|
1.000
|
12883.20
|
8691.54
|
|
Daily pressure experienced
|
0.95
|
0.89
|
1.01
|
1.001
|
17830.48
|
9885.92
|
|
Daily pressure utilized (partner’s view)
|
0.98
|
0.92
|
1.05
|
1.000
|
19278.32
|
10246.31
|
|
Daily pushing experienced
|
1.03
|
0.98
|
1.07
|
1.000
|
14025.59
|
9055.59
|
|
Daily pushing utilized (partner’s view)
|
1.03
|
0.99
|
1.07
|
1.000
|
21949.21
|
9623.34
|
|
Day
|
0.96
|
0.88
|
1.05
|
1.000
|
31435.94
|
8368.60
|
|
Daily weartime
|
1.00*
|
1.00
|
1.00
|
1.000
|
11933.20
|
7885.13
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
1.09
|
0.81
|
1.45
|
1.001
|
3608.13
|
6059.46
|
|
Mean persuasion utilized (partner’s view)
|
0.96
|
0.72
|
1.29
|
1.001
|
3610.97
|
5890.12
|
|
Mean pressure experienced
|
0.97
|
0.71
|
1.34
|
1.001
|
5511.94
|
7801.62
|
|
Mean pressure utilized (partner’s view)
|
0.98
|
0.72
|
1.32
|
1.001
|
4870.74
|
7220.18
|
|
Mean pushing experienced
|
1.02
|
0.66
|
1.56
|
1.000
|
5260.03
|
7662.48
|
|
Mean pushing utilized (partner’s view)
|
1.29
|
0.85
|
1.96
|
1.000
|
5271.31
|
7350.00
|
|
Mean weartime
|
1.00
|
1.00
|
1.00
|
1.000
|
17022.21
|
9961.31
|
|
Random Effects
|
|
sd(Intercept)
|
0.30
|
0.23
|
0.39
|
1.00
|
4214.82
|
7324.07
|
|
sd(Daily persuasion experienced)
|
0.05
|
0.03
|
0.08
|
1.00
|
7644.01
|
6901.55
|
|
sd(Daily persuasion utilized (partner’s view))
|
0.05
|
0.02
|
0.08
|
1.00
|
5351.19
|
3858.94
|
|
sd(Daily pressure experienced)
|
0.05
|
0.00
|
0.14
|
1.00
|
6659.89
|
6610.33
|
|
sd(Daily pressure utilized (partner’s view))
|
0.04
|
0.00
|
0.11
|
1.00
|
7701.49
|
6409.62
|
|
sd(Daily pushing experienced)
|
0.06
|
0.00
|
0.14
|
1.00
|
3294.00
|
5715.05
|
|
sd(Daily pushing utilized (partner’s view))
|
0.03
|
0.00
|
0.08
|
1.00
|
6464.76
|
7907.69
|
|
Additional Parameters
|
|
ar[1]
|
0.29
|
0.26
|
0.33
|
1.00
|
25732.94
|
7825.38
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sderr
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sigma
|
0.55
|
0.54
|
0.57
|
1.00
|
22829.74
|
8539.10
|
Affect
range(df_double$aff, na.rm = T)
## [1] 1 6
hist(df_double$aff, breaks = 15)

formula <- bf(
aff ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw | coupleID),
autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b"),
brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=6), # range of the outcome scale
brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),
brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = 5000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", "mood")
)









plot(pp_check(mood, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(mood, transform = log1p)
loo(pa_sub)
##
## Computed from 12000 by 3736 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -12060.6 177.6
## p_loo 30.2 2.8
## looic 24121.2 355.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.6]).
##
## Pareto k diagnostic values:
## Count Pct.
## (-Inf, 0.7] (good) 3732 99.9%
## (0.7, 1] (bad) 3 0.1%
## (1, Inf) (very bad) 1 0.0%
## Min. ESS
## (-Inf, 0.7] 662
## (0.7, 1] <NA>
## (1, Inf) <NA>
## See help('pareto-k-diagnostic') for details.
summarize_brms(
mood,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = F) %>%
print_df(rows_to_pack = rows_to_pack)
|
|
b
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
Intercept
|
4.73*
|
4.52
|
4.94
|
1.001
|
3046.66
|
4422.02
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
0.00
|
-0.03
|
0.04
|
1.000
|
17297.90
|
9842.88
|
|
Daily persuasion utilized (partner’s view)
|
0.02
|
-0.02
|
0.05
|
1.001
|
14451.84
|
9886.46
|
|
Daily pressure experienced
|
-0.05
|
-0.16
|
0.05
|
1.000
|
12202.80
|
8571.79
|
|
Daily pressure utilized (partner’s view)
|
-0.04
|
-0.18
|
0.09
|
1.001
|
11104.70
|
8944.62
|
|
Daily pushing experienced
|
0.02
|
-0.04
|
0.09
|
1.000
|
12680.12
|
9243.11
|
|
Daily pushing utilized (partner’s view)
|
0.06*
|
0.00
|
0.12
|
1.000
|
13534.43
|
9546.16
|
|
Day
|
0.22*
|
0.06
|
0.38
|
1.000
|
24416.16
|
8358.49
|
|
Daily weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
0.32
|
-0.26
|
0.88
|
1.001
|
2292.03
|
3556.47
|
|
Mean persuasion utilized (partner’s view)
|
0.23
|
-0.36
|
0.80
|
1.001
|
2249.94
|
3776.58
|
|
Mean pressure experienced
|
-0.29
|
-0.89
|
0.30
|
1.001
|
2937.39
|
4951.55
|
|
Mean pressure utilized (partner’s view)
|
-0.31
|
-0.89
|
0.28
|
1.001
|
2908.62
|
4194.29
|
|
Mean pushing experienced
|
0.25
|
-0.55
|
1.04
|
1.000
|
3993.12
|
6072.68
|
|
Mean pushing utilized (partner’s view)
|
0.39
|
-0.41
|
1.18
|
1.001
|
3968.95
|
6710.30
|
|
Mean weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Random Effects
|
|
sd(Intercept)
|
0.59
|
0.46
|
0.77
|
1.00
|
3439.67
|
6053.64
|
|
sd(Daily persuasion experienced)
|
0.03
|
0.00
|
0.07
|
1.00
|
5380.88
|
6051.41
|
|
sd(Daily persuasion utilized (partner’s view))
|
0.06
|
0.01
|
0.11
|
1.00
|
3288.91
|
3865.59
|
|
sd(Daily pressure experienced)
|
0.11
|
0.01
|
0.28
|
1.00
|
4275.51
|
5244.59
|
|
sd(Daily pressure utilized (partner’s view))
|
0.19
|
0.03
|
0.37
|
1.00
|
4126.88
|
4360.12
|
|
sd(Daily pushing experienced)
|
0.09
|
0.02
|
0.17
|
1.00
|
4973.05
|
3524.43
|
|
sd(Daily pushing utilized (partner’s view))
|
0.05
|
0.00
|
0.13
|
1.00
|
5294.70
|
5550.35
|
|
Additional Parameters
|
|
ar[1]
|
0.45
|
0.42
|
0.48
|
1.00
|
20444.09
|
8148.47
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sderr
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sigma
|
0.87
|
0.85
|
0.89
|
1.00
|
21206.93
|
8694.13
|
Reactance
Gaussian
range(df_double$reactance, na.rm = T)
## [1] 0 5
hist(df_double$reactance, breaks = 6)

formula <- bf(
reactance ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw | coupleID),
autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b"),
brms::set_prior("normal(0, 20)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),
brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1),
brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
reactance <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = 5000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", "reactance")
)
plot(reactance, ask = FALSE)









plot(pp_check(reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance, transform = log1p)
loo(reactance)
##
## Computed from 12000 by 756 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -1059.5 35.7
## p_loo 75.2 7.5
## looic 2118.9 71.3
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.8]).
##
## Pareto k diagnostic values:
## Count Pct.
## (-Inf, 0.7] (good) 747 98.8%
## (0.7, 1] (bad) 9 1.2%
## (1, Inf) (very bad) 0 0.0%
## Min. ESS
## (-Inf, 0.7] 174
## (0.7, 1] <NA>
## (1, Inf) <NA>
## See help('pareto-k-diagnostic') for details.
summarize_brms(
reactance,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = F) %>%
print_df(rows_to_pack = rows_to_pack)
|
|
b
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
Intercept
|
0.51*
|
0.33
|
0.70
|
1.000
|
16417.76
|
9892.94
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
-0.05
|
-0.11
|
0.01
|
1.000
|
19356.46
|
9444.65
|
|
Daily persuasion utilized (partner’s view)
|
0.00
|
-0.06
|
0.07
|
1.000
|
17176.01
|
9831.82
|
|
Daily pressure experienced
|
0.26*
|
0.04
|
0.47
|
1.000
|
11824.37
|
9336.31
|
|
Daily pressure utilized (partner’s view)
|
0.14
|
-0.06
|
0.40
|
1.000
|
9892.43
|
7761.10
|
|
Daily pushing experienced
|
0.08*
|
0.00
|
0.17
|
1.000
|
12338.51
|
9197.36
|
|
Daily pushing utilized (partner’s view)
|
-0.02
|
-0.10
|
0.07
|
1.000
|
19544.32
|
8481.36
|
|
Day
|
0.10
|
-0.17
|
0.35
|
1.000
|
24220.25
|
9140.35
|
|
Daily weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
0.08
|
-0.29
|
0.45
|
1.000
|
10206.81
|
9393.41
|
|
Mean persuasion utilized (partner’s view)
|
0.14
|
-0.25
|
0.55
|
1.000
|
10672.14
|
9527.96
|
|
Mean pressure experienced
|
0.62*
|
0.21
|
1.03
|
1.000
|
12564.34
|
9244.78
|
|
Mean pressure utilized (partner’s view)
|
0.16
|
-0.28
|
0.58
|
1.000
|
12489.27
|
9110.28
|
|
Mean pushing experienced
|
-0.28
|
-0.82
|
0.27
|
1.000
|
12347.15
|
10160.83
|
|
Mean pushing utilized (partner’s view)
|
-0.59*
|
-1.16
|
-0.01
|
1.000
|
13098.23
|
9577.68
|
|
Mean weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Random Effects
|
|
sd(Intercept)
|
0.27
|
0.15
|
0.41
|
1.00
|
6292.46
|
7778.91
|
|
sd(Daily persuasion experienced)
|
0.05
|
0.00
|
0.13
|
1.00
|
3879.53
|
6215.13
|
|
sd(Daily persuasion utilized (partner’s view))
|
0.04
|
0.00
|
0.13
|
1.00
|
7142.65
|
6976.51
|
|
sd(Daily pressure experienced)
|
0.39
|
0.22
|
0.62
|
1.00
|
6664.96
|
9032.62
|
|
sd(Daily pressure utilized (partner’s view))
|
0.25
|
0.01
|
0.63
|
1.00
|
3287.32
|
6408.13
|
|
sd(Daily pushing experienced)
|
0.10
|
0.01
|
0.23
|
1.00
|
4010.11
|
5667.82
|
|
sd(Daily pushing utilized (partner’s view))
|
0.05
|
0.00
|
0.15
|
1.00
|
8780.84
|
7897.48
|
|
Additional Parameters
|
|
ar[1]
|
0.01
|
-0.08
|
0.09
|
1.00
|
20581.26
|
9269.53
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sderr
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sigma
|
0.93
|
0.88
|
0.98
|
1.00
|
14809.72
|
8576.95
|
hypothesis(reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
## Hypothesis Estimate
## 1 (pressure_self_cw... > 0 0.18
## Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 0.12 -0.03 0.37 12.61
## Post.Prob Star
## 1 0.93
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Ordinal
df_double$reactance_ordinal <- factor(df_double$reactance,
levels = 0:5,
ordered = TRUE)
formula <- bf(
reactance_ordinal ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw | coupleID),
autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)
prior1 <- c(
set_prior("normal(0, 2.5)", class = "b"),
set_prior("normal(0, 5)", class = "Intercept"),
set_prior("normal(0, 1)", class = "sd", group = "coupleID", lb = 0),
set_prior("normal(0, 0.075)", class = "ar", lb = -1, ub = 1),
set_prior("normal(0.5, 2.0)", class = "sderr", lb = 0)
)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
reactance_ordinal <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::cumulative(),
control = list(adapt_delta = 0.95),
iter = 5000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", "reactance_ordinal")
)
plot(reactance_ordinal, ask = FALSE)










plot(pp_check(reactance_ordinal, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(reactance_ordinal))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(reactance_ordinal, transform = log1p)
loo(reactance_ordinal)
##
## Computed from 12000 by 756 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -683.1 32.1
## p_loo 112.0 7.4
## looic 1366.2 64.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.1, 1.3]).
##
## Pareto k diagnostic values:
## Count Pct.
## (-Inf, 0.7] (good) 753 99.6%
## (0.7, 1] (bad) 3 0.4%
## (1, Inf) (very bad) 0 0.0%
## Min. ESS
## (-Inf, 0.7] 143
## (0.7, 1] <NA>
## (1, Inf) <NA>
## See help('pareto-k-diagnostic') for details.
summarize_brms(
reactance_ordinal,
model_rows_fixed = model_rows_fixed_ordinal,
model_rows_random = model_rows_random_ordinal,
model_rownames_fixed = model_rownames_fixed_ordinal,
model_rownames_random = model_rownames_random_ordinal,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack_ordinal)
|
|
OR
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
Intercepts
|
|
Intercept
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Intercept[1]
|
4.29*
|
2.49
|
8.21
|
1.001
|
2256.01
|
1951.78
|
|
Intercept[2]
|
9.75*
|
5.40
|
21.51
|
1.002
|
1435.82
|
1081.80
|
|
Intercept[3]
|
28.61*
|
14.76
|
75.24
|
1.003
|
1210.60
|
778.30
|
|
Intercept[4]
|
133.42*
|
58.72
|
466.56
|
1.004
|
1094.69
|
724.17
|
|
Intercept[5]
|
5469.93*
|
1336.97
|
45257.45
|
1.003
|
1166.93
|
814.31
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
0.84*
|
0.70
|
0.99
|
1.000
|
6298.44
|
5973.36
|
|
Daily persuasion utilized (partner’s view)
|
1.02
|
0.83
|
1.24
|
1.000
|
10312.32
|
8638.15
|
|
Daily pressure experienced
|
1.92*
|
1.20
|
2.93
|
1.001
|
3247.00
|
2869.57
|
|
Daily pressure utilized (partner’s view)
|
1.25
|
0.74
|
2.11
|
1.000
|
8335.13
|
7024.63
|
|
Daily pushing experienced
|
1.17
|
0.96
|
1.46
|
1.000
|
7133.18
|
6701.14
|
|
Daily pushing utilized (partner’s view)
|
0.91
|
0.69
|
1.19
|
1.001
|
9296.13
|
7955.74
|
|
Day
|
1.48
|
0.72
|
3.06
|
1.000
|
11603.05
|
7629.89
|
|
Daily weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
1.11
|
0.38
|
3.29
|
1.001
|
5666.63
|
6811.76
|
|
Mean persuasion utilized (partner’s view)
|
1.40
|
0.45
|
4.79
|
1.001
|
5800.97
|
6413.49
|
|
Mean pressure experienced
|
3.71*
|
1.20
|
12.26
|
1.001
|
5312.91
|
5602.60
|
|
Mean pressure utilized (partner’s view)
|
1.20
|
0.36
|
3.84
|
1.001
|
6106.53
|
7633.92
|
|
Mean pushing experienced
|
1.18
|
0.26
|
5.67
|
1.000
|
7661.62
|
8475.36
|
|
Mean pushing utilized (partner’s view)
|
0.09*
|
0.01
|
0.57
|
1.001
|
5985.19
|
7322.95
|
|
Mean weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Random Effects
|
|
sd(Intercept)
|
0.84
|
0.48
|
1.31
|
1.00
|
2922.60
|
3336.30
|
|
sd(Daily persuasion experienced)
|
0.18
|
0.01
|
0.44
|
1.00
|
2882.76
|
5395.24
|
|
sd(Daily persuasion utilized (partner’s view))
|
0.22
|
0.01
|
0.53
|
1.00
|
3661.25
|
5972.87
|
|
sd(Daily pressure experienced)
|
0.57
|
0.08
|
1.15
|
1.00
|
2764.88
|
2568.59
|
|
sd(Daily pressure utilized (partner’s view))
|
0.45
|
0.02
|
1.34
|
1.00
|
4034.58
|
7034.00
|
|
sd(Daily pushing experienced)
|
0.22
|
0.01
|
0.52
|
1.00
|
3442.46
|
5129.78
|
|
sd(Daily pushing utilized (partner’s view))
|
0.19
|
0.01
|
0.61
|
1.00
|
4713.67
|
6572.87
|
|
Additional Parameters
|
|
ar[1]
|
0.00
|
-0.15
|
0.14
|
1.00
|
16982.67
|
8940.39
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sderr
|
0.51
|
0.02
|
1.65
|
1.01
|
613.02
|
433.90
|
|
sigma
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
disc
|
1.00
|
1.00
|
1.00
|
NA
|
NA
|
NA
|
hypothesis(reactance_ordinal, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
## Hypothesis Estimate
## 1 (pressure_self_cw... > 0 0.49
## Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 0.26 0.05 0.89 26.84
## Post.Prob Star
## 1 0.96 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Binary
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}
formula <- bf(
is_reactance ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
pushing_self_cb + pushing_partner_cb +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
pushing_self_cw + pushing_partner_cw | coupleID),
autocor = ~ ar(time = day, gr = coupleID:userID, p = 1)
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b"),
brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5), # range of the outcome scale
brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0),
brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1)
#brms::set_prior("cauchy(0, 10)", class = "sigma", lb = 0)
)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = 5000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", "is_reactance")
)
plot(is_reactance, ask = FALSE)









plot(pp_check(is_reactance, type = 'ecdf_overlay'))
## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

plot(pp_check(is_reactance))
## Using 10 posterior draws for ppc type 'dens_overlay' by default.

#pp_check_transformed(is_reactance, transform = log1p)
loo(is_reactance)
##
## Computed from 12000 by 756 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -328.4 12.7
## p_loo 260.2 10.7
## looic 656.8 25.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.2]).
##
## Pareto k diagnostic values:
## Count Pct.
## (-Inf, 0.7] (good) 35 4.6%
## (0.7, 1] (bad) 572 75.7%
## (1, Inf) (very bad) 149 19.7%
## Min. ESS
## (-Inf, 0.7] 396
## (0.7, 1] <NA>
## (1, Inf) <NA>
## See help('pareto-k-diagnostic') for details.
summarize_brms(
is_reactance,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)
|
|
OR
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
Intercept
|
0.03*
|
0.00
|
0.18
|
1.001
|
2961.20
|
5776.43
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
0.55
|
0.26
|
1.01
|
1.000
|
4130.22
|
6033.87
|
|
Daily persuasion utilized (partner’s view)
|
1.86
|
0.76
|
6.03
|
1.001
|
3737.98
|
5220.40
|
|
Daily pressure experienced
|
8.70*
|
1.69
|
60.77
|
1.000
|
3710.48
|
6299.81
|
|
Daily pressure utilized (partner’s view)
|
2.19
|
0.33
|
18.42
|
1.000
|
6729.25
|
6651.90
|
|
Daily pushing experienced
|
2.34*
|
1.07
|
6.63
|
1.000
|
3633.87
|
4827.24
|
|
Daily pushing utilized (partner’s view)
|
0.72
|
0.21
|
2.27
|
1.001
|
6444.05
|
6815.31
|
|
Day
|
3.78
|
0.38
|
43.96
|
1.000
|
6637.82
|
7554.77
|
|
Daily weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
13.42
|
0.90
|
291.21
|
1.000
|
3985.43
|
6196.30
|
|
Mean persuasion utilized (partner’s view)
|
5.05
|
0.31
|
91.34
|
1.000
|
7444.78
|
7562.53
|
|
Mean pressure experienced
|
132.06*
|
3.62
|
5363.94
|
1.000
|
9338.38
|
8352.40
|
|
Mean pressure utilized (partner’s view)
|
6.03
|
0.15
|
270.01
|
1.000
|
8879.00
|
9243.83
|
|
Mean pushing experienced
|
4.47
|
0.10
|
220.19
|
1.000
|
7718.10
|
8220.27
|
|
Mean pushing utilized (partner’s view)
|
0.06
|
0.00
|
3.21
|
1.000
|
8490.07
|
8668.77
|
|
Mean weartime
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Random Effects
|
|
sd(Intercept)
|
3.39
|
1.86
|
5.30
|
1.00
|
2132.54
|
3149.41
|
|
sd(Daily persuasion experienced)
|
0.58
|
0.03
|
1.50
|
1.00
|
1978.49
|
4478.41
|
|
sd(Daily persuasion utilized (partner’s view))
|
1.50
|
0.42
|
2.95
|
1.00
|
2329.46
|
2696.99
|
|
sd(Daily pressure experienced)
|
1.97
|
0.11
|
4.47
|
1.00
|
1972.68
|
3600.22
|
|
sd(Daily pressure utilized (partner’s view))
|
1.35
|
0.05
|
3.81
|
1.00
|
4817.26
|
5523.96
|
|
sd(Daily pushing experienced)
|
0.84
|
0.05
|
2.02
|
1.00
|
2576.31
|
3993.09
|
|
sd(Daily pushing utilized (partner’s view))
|
0.75
|
0.03
|
2.23
|
1.00
|
4299.27
|
6263.29
|
|
Additional Parameters
|
|
ar[1]
|
0.06
|
-0.16
|
0.27
|
1.00
|
1980.35
|
3211.48
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sderr
|
5.90
|
3.10
|
9.31
|
1.00
|
1442.76
|
1679.03
|
|
sigma
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
## Hypothesis Estimate
## 1 (pressure_self_cw... > 0 1.31
## Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 0.99 -0.23 3.01 11.11
## Post.Prob Star
## 1 0.92
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
Report All Models
if (report_ordinal) {
model_rows_random_final <- model_rows_random_ordinal
model_rows_fixed_final <- model_rows_fixed_ordinal
model_rownames_fixed_final <- model_rownames_fixed_ordinal
model_rownames_random_final <- model_rownames_random_ordinal
rows_to_pack_final <- rows_to_pack_ordinal
} else {
model_rows_random_final <- model_rows_random
model_rows_fixed_final <- model_rows_fixed
model_rownames_fixed_final <- model_rownames_fixed
model_rownames_random_final <- model_rownames_random
rows_to_pack_final <- rows_to_pack
}
all_models <- report_side_by_side(
pa_sub,
pa_obj_log,
mood,
reactance,
is_reactance,
model_rows_random = model_rows_random_final,
model_rows_fixed = model_rows_fixed_final,
model_rownames_random = model_rownames_random_final,
model_rownames_fixed = model_rownames_fixed_final
)
## [1] "pa_sub"
## [1] "pa_obj_log"
## [1] "mood"
## [1] "reactance"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
print_df(rows_to_pack = rows_to_pack_final) %>%
add_header_above(
c(" ", "Subjective MVPA" = 2,
"Device-Based MVPA" = 2,
"Mood" = 2,
"Reactance Gaussian" = 2,
"Reactance Dichotome" = 2)
)
export_xlsx(
summary_all_models,
rows_to_pack = rows_to_pack_final,
file.path("Output", "AllModels.xlsx"),
merge_option = 'both',
simplify_2nd_row = TRUE,
colwidths = c(38, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
line_above_rows = c(1,2),
line_below_rows = c(-1)
)
summary_all_models
|
|
Subjective MVPA
|
Device-Based MVPA
|
Mood
|
Reactance Gaussian
|
Reactance Dichotome
|
|
|
IRR pa_sub
|
95% CI pa_sub
|
exp(Est.) pa_obj_log
|
95% CI pa_obj_log
|
b mood
|
95% CI mood
|
b reactance
|
95% CI reactance
|
OR is_reactance
|
95% CI is_reactance
|
|
Intercept
|
27.35*
|
[21.03, 35.95]
|
118.09*
|
[106.05, 131.83]
|
4.73*
|
[ 4.52, 4.94]
|
0.51*
|
[ 0.33, 0.70]
|
0.03*
|
[0.00, 0.18]
|
|
Within-Person Effects
|
|
Daily persuasion experienced
|
1.20*
|
[ 1.07, 1.35]
|
1.03
|
[ 1.00, 1.06]
|
0.00
|
[-0.03, 0.04]
|
-0.05
|
[-0.11, 0.01]
|
0.55
|
[0.26, 1.01]
|
|
Daily persuasion utilized (partner’s view)
|
1.19*
|
[ 1.06, 1.34]
|
1.02
|
[ 0.99, 1.05]
|
0.02
|
[-0.02, 0.05]
|
0.00
|
[-0.06, 0.07]
|
1.86
|
[0.76, 6.03]
|
|
Daily pressure experienced
|
0.94
|
[ 0.71, 1.30]
|
0.95
|
[ 0.89, 1.01]
|
-0.05
|
[-0.16, 0.05]
|
0.26*
|
[ 0.04, 0.47]
|
8.70*
|
[1.69, 60.77]
|
|
Daily pressure utilized (partner’s view)
|
1.16
|
[ 0.88, 1.59]
|
0.98
|
[ 0.92, 1.05]
|
-0.04
|
[-0.18, 0.09]
|
0.14
|
[-0.06, 0.40]
|
2.19
|
[0.33, 18.42]
|
|
Daily pushing experienced
|
1.13
|
[ 0.92, 1.41]
|
1.03
|
[ 0.98, 1.07]
|
0.02
|
[-0.04, 0.09]
|
0.08*
|
[ 0.00, 0.17]
|
2.34*
|
[1.07, 6.63]
|
|
Daily pushing utilized (partner’s view)
|
1.13
|
[ 0.95, 1.36]
|
1.03
|
[ 0.99, 1.07]
|
0.06*
|
[ 0.00, 0.12]
|
-0.02
|
[-0.10, 0.07]
|
0.72
|
[0.21, 2.27]
|
|
Day
|
0.79
|
[ 0.58, 1.10]
|
0.96
|
[ 0.88, 1.05]
|
0.22*
|
[ 0.06, 0.38]
|
0.10
|
[-0.17, 0.35]
|
3.78
|
[0.38, 43.96]
|
|
Daily weartime
|
NA
|
NA
|
1.00*
|
[ 1.00, 1.00]
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Between-Person Effects
|
|
Mean persuasion experienced
|
1.61
|
[ 0.83, 3.13]
|
1.09
|
[ 0.81, 1.45]
|
0.32
|
[-0.26, 0.88]
|
0.08
|
[-0.29, 0.45]
|
13.42
|
[0.90, 291.21]
|
|
Mean persuasion utilized (partner’s view)
|
1.19
|
[ 0.62, 2.28]
|
0.96
|
[ 0.72, 1.29]
|
0.23
|
[-0.36, 0.80]
|
0.14
|
[-0.25, 0.55]
|
5.05
|
[0.31, 91.34]
|
|
Mean pressure experienced
|
0.50
|
[ 0.22, 1.11]
|
0.97
|
[ 0.71, 1.34]
|
-0.29
|
[-0.89, 0.30]
|
0.62*
|
[ 0.21, 1.03]
|
132.06*
|
[3.62, 5363.94]
|
|
Mean pressure utilized (partner’s view)
|
0.43*
|
[ 0.19, 0.96]
|
0.98
|
[ 0.72, 1.32]
|
-0.31
|
[-0.89, 0.28]
|
0.16
|
[-0.28, 0.58]
|
6.03
|
[0.15, 270.01]
|
|
Mean pushing experienced
|
1.71
|
[ 0.64, 4.56]
|
1.02
|
[ 0.66, 1.56]
|
0.25
|
[-0.55, 1.04]
|
-0.28
|
[-0.82, 0.27]
|
4.47
|
[0.10, 220.19]
|
|
Mean pushing utilized (partner’s view)
|
2.03
|
[ 0.76, 5.67]
|
1.29
|
[ 0.85, 1.96]
|
0.39
|
[-0.41, 1.18]
|
-0.59*
|
[-1.16, -0.01]
|
0.06
|
[0.00, 3.21]
|
|
Mean weartime
|
NA
|
NA
|
1.00
|
[ 1.00, 1.00]
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
Random Effects
|
|
sd(Intercept)
|
0.61
|
[ 0.45, 0.83]
|
0.30
|
[0.23, 0.39]
|
0.59
|
[0.46, 0.77]
|
0.27
|
[ 0.15, 0.41]
|
3.39
|
[ 1.86, 5.30]
|
|
sd(Daily persuasion experienced)
|
0.09
|
[ 0.00, 0.24]
|
0.05
|
[0.03, 0.08]
|
0.03
|
[0.00, 0.07]
|
0.05
|
[ 0.00, 0.13]
|
0.58
|
[ 0.03, 1.50]
|
|
sd(Daily persuasion utilized (partner’s view))
|
0.08
|
[ 0.00, 0.21]
|
0.05
|
[0.02, 0.08]
|
0.06
|
[0.01, 0.11]
|
0.04
|
[ 0.00, 0.13]
|
1.50
|
[ 0.42, 2.95]
|
|
sd(Daily pressure experienced)
|
0.17
|
[ 0.01, 0.53]
|
0.05
|
[0.00, 0.14]
|
0.11
|
[0.01, 0.28]
|
0.39
|
[ 0.22, 0.62]
|
1.97
|
[ 0.11, 4.47]
|
|
sd(Daily pressure utilized (partner’s view))
|
0.16
|
[ 0.01, 0.49]
|
0.04
|
[0.00, 0.11]
|
0.19
|
[0.03, 0.37]
|
0.25
|
[ 0.01, 0.63]
|
1.35
|
[ 0.05, 3.81]
|
|
sd(Daily pushing experienced)
|
0.28
|
[ 0.02, 0.58]
|
0.06
|
[0.00, 0.14]
|
0.09
|
[0.02, 0.17]
|
0.10
|
[ 0.01, 0.23]
|
0.84
|
[ 0.05, 2.02]
|
|
sd(Daily pushing utilized (partner’s view))
|
0.12
|
[ 0.00, 0.32]
|
0.03
|
[0.00, 0.08]
|
0.05
|
[0.00, 0.13]
|
0.05
|
[ 0.00, 0.15]
|
0.75
|
[ 0.03, 2.23]
|
|
Additional Parameters
|
|
ar[1]
|
0.03
|
[-0.94, 0.94]
|
0.29
|
[0.26, 0.33]
|
0.45
|
[0.42, 0.48]
|
0.01
|
[-0.08, 0.09]
|
0.06
|
[-0.16, 0.27]
|
|
nu
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
shape
|
0.14
|
[ 0.13, 0.14]
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
|
sderr
|
0.05
|
[ 0.00, 0.13]
|
NA
|
NA
|
NA
|
NA
|
NA
|
NA
|
5.90
|
[ 3.10, 9.31]
|
|
sigma
|
NA
|
NA
|
0.55
|
[0.54, 0.57]
|
0.87
|
[0.85, 0.89]
|
0.93
|
[ 0.88, 0.98]
|
NA
|
NA
|
Analyses were conducted using the R Statistical language (version
4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)
- Allaire J, Xie Y, Dervieux C, McPherson J, Luraschi J, Ushey K,
Atkins A, Wickham H, Cheng J, Chang W, Iannone R (2024). rmarkdown:
Dynamic Documents for R. R package version 2.28, https://github.com/rstudio/rmarkdown. Xie Y, Allaire J,
Grolemund G (2018). R Markdown: The Definitive Guide. Chapman
and Hall/CRC, Boca Raton, Florida. ISBN 9781138359338, https://bookdown.org/yihui/rmarkdown. Xie Y, Dervieux C,
Riederer E (2020). R Markdown Cookbook. Chapman and Hall/CRC,
Boca Raton, Florida. ISBN 9780367563837, https://bookdown.org/yihui/rmarkdown-cookbook.
- Bååth R (2024). beepr: Easily Play Notification Sounds on any
Platform. R package version 2.0, https://CRAN.R-project.org/package=beepr.
- Bengtsson H (2003). “The R.oo package - Object-Oriented Programming
with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A
(eds.), Proceedings of the 3rd International Workshop on Distributed
Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
- Bengtsson H (2003). “The R.oo package - Object-Oriented Programming
with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A
(eds.), Proceedings of the 3rd International Workshop on Distributed
Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
- Bengtsson H (2023). R.utils: Various Programming Utilities.
R package version 2.12.3, https://CRAN.R-project.org/package=R.utils.
- Bürkner P (2017). “brms: An R Package for Bayesian Multilevel Models
Using Stan.” Journal of Statistical Software, 80(1),
1-28. doi:10.18637/jss.v080.i01 https://doi.org/10.18637/jss.v080.i01. Bürkner P (2018).
“Advanced Bayesian Multilevel Modeling with the R Package brms.” The
R Journal, 10(1), 395-411. doi:10.32614/RJ-2018-017
https://doi.org/10.32614/RJ-2018-017. Bürkner P (2021).
“Bayesian Item Response Modeling in R with brms and Stan.” Journal
of Statistical Software, 100(5), 1-54. doi:10.18637/jss.v100.i05 https://doi.org/10.18637/jss.v100.i05.
- Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N,
Ucar I, Bates D, Chambers J (2024). Rcpp: Seamless R and C++
Integration. R package version 1.0.13, https://CRAN.R-project.org/package=Rcpp. Eddelbuettel D,
François R (2011). “Rcpp: Seamless R and C++ Integration.” Journal
of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08 https://doi.org/10.18637/jss.v040.i08. Eddelbuettel D
(2013). Seamless R and C++ Integration with Rcpp. Springer, New
York. doi:10.1007/978-1-4614-6868-4 https://doi.org/10.1007/978-1-4614-6868-4, ISBN
978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). “Extending R with
C++: A Brief Introduction to Rcpp.” The American Statistician,
72(1), 28-36. doi:10.1080/00031305.2017.1375990 https://doi.org/10.1080/00031305.2017.1375990.
- Gabry J, Mahr T (2024). “bayesplot: Plotting for Bayesian Models.” R
package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry J, Simpson D,
Vehtari A, Betancourt M, Gelman A (2019). “Visualization in Bayesian
workflow.” J. R. Stat. Soc. A, 182, 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378.
- Grolemund G, Wickham H (2011). “Dates and Times Made Easy with
lubridate.” Journal of Statistical Software, 40(3),
1-25. https://www.jstatsoft.org/v40/i03/.
- Hartig F (2022). DHARMa: Residual Diagnostics for Hierarchical
(Multi-Level / Mixed) Regression Models. R package version 0.4.6,
https://CRAN.R-project.org/package=DHARMa.
- Hester J, Wickham H, Csárdi G (2024). fs: Cross-Platform File
System Operations Based on ‘libuv’. R package version 1.6.4, https://CRAN.R-project.org/package=fs.
- Küng P (2023). wbCorr: Bivariate Within- and Between-Cluster
Correlations. University of Zürich. R package version 0.1.22. https://github.com/Pascal-Kueng/wbCorr.
- Müller K, Wickham H (2023). tibble: Simple Data Frames. R
package version 3.2.1, https://CRAN.R-project.org/package=tibble.
- R Core Team (2024). R: A Language and Environment for
Statistical Computing. R Foundation for Statistical Computing,
Vienna, Austria. https://www.R-project.org/.
- Schauberger P, Walker A (2024). openxlsx: Read, Write and Edit
xlsx Files. R package version 4.2.7.1, https://CRAN.R-project.org/package=openxlsx.
- Wickham H (2016). ggplot2: Elegant Graphics for Data
Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
- Wickham H (2023). forcats: Tools for Working with Categorical
Variables (Factors). R package version 1.0.0, https://CRAN.R-project.org/package=forcats.
- Wickham H (2023). stringr: Simple, Consistent Wrappers for
Common String Operations. R package version 1.5.1, https://CRAN.R-project.org/package=stringr.
- Wickham H (2024). rvest: Easily Harvest (Scrape) Web Pages.
R package version 1.0.4, https://CRAN.R-project.org/package=rvest.
- Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E,
Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K,
Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.”
Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686
https://doi.org/10.21105/joss.01686.
- Wickham H, Bryan J (2023). readxl: Read Excel Files. R
package version 1.4.3, https://CRAN.R-project.org/package=readxl.
- Wickham H, François R, Henry L, Müller K, Vaughan D (2023).
dplyr: A Grammar of Data Manipulation. R package version 1.1.4,
https://CRAN.R-project.org/package=dplyr.
- Wickham H, Henry L (2023). purrr: Functional Programming
Tools. R package version 1.0.2, https://CRAN.R-project.org/package=purrr.
- Wickham H, Hester J, Bryan J (2024). readr: Read Rectangular
Text Data. R package version 2.1.5, https://CRAN.R-project.org/package=readr.
- Wickham H, Hester J, Ooms J (2023). xml2: Parse XML. R
package version 1.3.6, https://CRAN.R-project.org/package=xml2.
- Wickham H, Vaughan D, Girlich M (2024). tidyr: Tidy Messy
Data. R package version 1.3.1, https://CRAN.R-project.org/package=tidyr.
- Xie Y (2024). knitr: A General-Purpose Package for Dynamic
Report Generation in R. R package version 1.48, https://yihui.org/knitr/. Xie Y (2015). Dynamic
Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca
Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A
Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch
F, Peng RD (eds.), Implementing Reproducible Computational
Research. Chapman and Hall/CRC. ISBN 978-1466561595.
- Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’
and Pipe Syntax. R package version 1.4.0, https://CRAN.R-project.org/package=kableExtra.